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SECTION  I 


INTRODUCTION 


MICE  is  an  implicit,  time-step-split,  MHD  code  which  is  being 
developed  at  Mission  Research  Corporation  to  extend  and  improve  upon 
our  capability  to  compute  the  late-time  phenomenology  of  high  altitude 
nuclear  explosions.  The  code  has  been  designed  to  do  calculations  in 
one-dimensional  spherical,  1-D  planar,  2-D  cartesian,  2-D  cylindrical, 
or  3-D  cartesian  geometry. 

In  addition  to  computing  an  approximate  numerical  solution  to 
the  one-fluid  MHD  equations,  the  code  also  calculates  the  following 
for  each  space  point  and  for  each  time,  (1)  the  Lagrangian  coordinate, 

(2)  the  number  densities  of  the  atmospheric  species  N2>  07,  N,  0,  N+, 

0+,  and  an  "all-purpose"  molecular  ion  X*  ,  and  (3)  the  rate  of  energy 
loss  by  radiation  of  UV.  There  is  a  deposition  routine  which  may  be 
called  one  or  more  times  during  the  running  of  a  problem,  to  deposit 
x-rays  and/or  UV  radiation.  A  very  general  rezoning  routine  allows  a 
problem  to  be  transferred  to  any  new  grid  at  any  time  during  the  run¬ 
ning  of  the  problem.  These  "nW^nomenology"  aspects  of  the  code  have 
been  patterned  after  the  MP.C  code  HOIL,  which  has  been  described  in 
reference  1 . 

The  DNA  sponsored  high-altitudo  phenomenology  community  already 
has  several  codes  which  calculate!  late-time  phenomenology.  There  are 
two-dimensional  codes  at  LASL,  NRL,  and  MRC,  and  three-dimensional  codes 
at  AFWL  and  NRL.  These  codes  all  use  explicit  difference  techniques. 

The  characteristic  feature  of  such  difference  techniques  is  that 
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numerical  stability  considerations  impose  the  condition  that  the  mesh 
speed,  Ax/At,  must  exceed  the  physical  signal  propagation  speed,  U  +  C, 
where  U  is  the  magnitude  of  the  fluid  velocity,  and  C  is  the  speed 
with  which  waves  travel  through  the  fluid.  The  corresponding  require¬ 
ment  for  implicit  codes,  such  as  MICE,  is  simply  Ax/At  >  U  .  For 
strong  shock  problems,  where  U  and  C  are  comparable  in  magnitude, 
the  two  type  of  codes  require  comparable  running  times.  However,  for 
problems  involving  subsonic  flow,  such  as  low  altitude  fireball  rise, 
or  the  late-time  stages  ci  high  altitude  atmospheric  heave,  the  implicit 
type  of  code  will  be  able  to  use  much  larger  time  steps,  and  therefore, 
require  much  less  computer  time.  It  is  the  need  for  shorten  running 
times  for  high  altitude  heave  calculations  which  has  prompted  the 
development  of  MICE,  In  the  ROSC  code  (Radar  and  Optica]  Systems  Code) 
it  will,  of  course,  be  necessary  to  have  a  reliable  model  of  atmospheric 
heave.  One  possibility  for  this  "model",  is  to  use  a  coarsely  zoned 
3D  MHD  code  as  an  integral  part  of  ROSC.  If  this  is  to  be  done,  the 
MHD  scheme  which  is  used  will  need  to  be  as  efficient  and  fast  running 
as  possible.  An  implicit  code,  such  as  MICE,  is  needed  to  meet  these 
requirements . 

Another  distinctive  feature  of  MICE  is  its  use  of  the 
approximation  of  time-step-splitting.  This  greatly  simplifies  the 
coding,  makes  it  easier  to  implement  changes  in  the  difference  scheme, 
reduces  the  main  core  requirements,  and  gives  the  code  the  ability  to 
do  calculations  in  several  different  geometries.  This  aspect  of  the 
code  is  discussed  further  in  section  3  of  this  report,  which  deals  with 
the  details  of  the  difference  scheme. 
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SECTION  II 


SOME  NUMERICAL  RESULTS 

Here  we  present  the  results  of  two  test  calculations.  The 
first  calculation  is  that  of  the  early  expansion  and  rise  of  a  hot 
bubble  in  an  exponential  atmosphere.  Some  comparisons  are  made  with 
other  numerical  solutions  of  the  same  problem.  The  second  calculation 
is  that  of  the  buoyant  rise  o^  a  warm  cylinder. 

HOT  BUBBLE  PROBLEM 

This  problem  was  chosen  as  a  test  case,  because  numerical  solutions 

(21 

of  the  same  problem  have  been  published  by  Harlow  and  Meixne1'  and 
by  Sappenfield^ .  The  initial  conditions  are  as  follows:  8.2  x  102°  ergs 
of  energy  are  placed  within  a  sphere  of  2.7  km  radius.  The  velocity 
within  the  sphere  is  taken  to  be  directed  radially  outward  from  the 
center,  and  to  have  the  magnitude  v  ■  4.25  1  +  j  km/sec,  in  which 

Ah  is  the  vertical  displacement  from  the  burst  point.  The  internal 
energy  within  the  sphere  is  8  x  1012  ergs/gm.  Burst  point  mass  density 
is  1.2  x  10  s  gm/cc,  and  the  atmospheric,  scale  height  is  6  km.  In  the 
MICE  calculation,  the  grid  spacing  was  initially  chosen  as  AR  ■  AZ  «=  .4  km. 

At  t  »  .45  sec,  the  problem  was  resoned  into  a  grid  with  AR  «*  AZ  ■  1  km. 

A  numerical  solution  to  this  problem  was  obtained  by  Harlow 
and  Meixner  with  a  particle-in-cell  (PIC)  code,  and  later  by  Sappenfield 
with  a  mixed  Eulerian-lagrangian  code  called  H2SE.  Figure  1  shows  the 
shock  position  as  a  function  of  time  for  the  three  directions,  "down", 
"sideways",  and  "up",  as  calculated  by  PIC,  H2SE,  and  MICE.  Of  the 
three  calculations,  the  H2SE  result  is  probably  the  most  reliable, 


TIME  (SEC) 

Figure  1.  Shock  positions  as  a  function 
of  time  for  the  hot  bubble 
problem. 
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since  H2SE  has  Lagrangian  boundaries  perpendicular  to  the  direction  of 
shock  propagation.  Lagrangian  codes  usually  give  better  definition  of 
shock  positions  than  do  Eulerian  codes. 

Figure  2  shows  contour  of  mass  density  as  given  by  H2SE  and 
MICE.  The  two  sets  of  contours  agree  qualitatively  very  well.  The 
lack  of  detailed  agreement  is  caused  by  the  fact  that  MICE  is  an  Eulerian 
code,  and  used  a  mesh  spacing  of  1  km  for  this  problem. 

BUOYANT  RISE  OF  A  WARM  CYLINDER 

Recently,  at  MRC,  an  analytic  solution  has  been  obtained  for 

the  early  buoyant  rise  of  a  cylinder  of  imcompressible  fluid  of  mass 

density  pj  ,  imbedded  in  another  fluid  of  density  p0  >  pi  ,  in  a 
4 

gravitational  field  .  This  solution  was  obtained  by  expanding  the  fluid 
equations  for  inviscid,  imcompressible  flow,  in  a  power  series  in  time, 
and  is  an  exact  solution  for  as  long  as  the  power  series  converges. 

We  give  here  a  comparison  of  the  analytic  solution  for  the  case  p*  *  .$p0, 
and  a  similar  problem  which  was  run  with  the  MICE  code. 

The  initial  conditions  for  the  MICE  problem  were  chosen  as 

fellows : 

ambient  mass  density  p0  B  8.3x10  5  gm/cc 

density  scale  height  Hg  ■  6.3xl05  cm 

ambient  specific  internal  energy  I 5  «  1.53x10®  org/gm 

spocific  heat  ratio  y  ■  1.4 

initial  cylinder  radius  R  *  5xl04  cm 

Inside  the  cylinder  the  mass  density  wa'  tken  to  be  .5  x  ambient  density, 
and  the  internal  energv  2.0  x  ambient  internal  energy.  The  problem  was 
run  in  2-D  Cartesian  geometry,  with  a  minimum  cell  size  of  about  .1  R. 
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lass  density  at  t  =  2.85  sec  for  the  hot  bubble  problem 


t 


This  problem,  as  done  by  the  MICE  code,  is,  of  course,  a 
compressible  flow  problem,  since  the  sound  speed  is  finite.  Also,  since 
this  code,  like  all  fluid  codes,  has  inherent  numerical  viscous  effects, 
the  problem  being  reported  here  is  clearly  a  different  problem  than 
that  for  which  we  have  the  analytic  solution.  However,  if  the  sound 
speed  is  sufficiently  large,  and  if  the  numerical  viscous  effects  do 
not  dominate  the  problem,  we  would  expect  the  results  to  be  not  too 
different  from  the  analytic  results.  To  see  how  large  the  sound  speed 
is  in  the  context  of  this  problem,  it  is  useful  to  work  in  an  appropriate 
set  of  dimensionless  units. 


The  power  series  expansion  for  the  height  of  the  center  of 
mass  of  the  cylinder  as  a  function  of  time  is 


h  *  . r  g  t!  ♦  higher  order  terms, 

Pc  +  Pi  <■ 


which  may  be  written  in  dimensionless  form  as 


H  - 


h 

—  a 

R 


V 


T2  + 


where  g' 


ftA.-JU.  r 

Po  ♦  Pi 


R  is  the  initial  radius  of  the  cylinder,  and  g  is  the  acceleration  of 
gravity.  If  c  is  the  sound  speed  in  cgs  units,  then  the  dimensionless 
sound  speed  is: 

C  *  c  x  ~  x  v/R/g '  ■  c/\~R  U 1 


For  the  set  of  initial  conditions  given  above,  the  sound  speed  inside 
the  cylinder  is 


C 


[y(Y-l)x2I0/(R 


10.3 
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Thus,  a  sound  signal  can  traverse  the  radius  of  the  cylinder  about  10 
times  in  the  time  required  fur  the  cylinder  to  rise  a  distance  of  i-  R  . 
We  would  therefore  expect  "compressibility"  effects  to  be  relatively 
minor.  The  numerical  viscosity  effects  are  more  difficult  to  assess, 
and  we  have  made  no  attempts  to  do  so. 

Figure  3  shows  a  comparison  of  the  position  of  the  boundary 
of  the  cylinder  at  T2  =  1.5.  The  analytic  solution  fails  to  converge 
for  times  much  beyond  this  time.  The  MICE  boundary  agrees  reasonably 
well  with  the  analytic  solution,  although  it  tends  to  oscillate  about 
the  analytic  result.  Figure  4  shows  the  best  MICE  calculation  we  could 
get  for  this  problem.  Here,  we  used  a  minimum  cell  spacing  of  about 
.05  R,  with  a  corresponding  decrease  in  At. 
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MICE  Calculation  of  the  Buoyant  Rise  of 
a  Warm  Cylinder,  with  Minimum  AX  R/10. 


Altitude 


Figure  4.  MICE  Calculation  of  the  Buoyant  Rise  of 
a  Warm  Cylinder,  with  Minimum  AX  »  F/20 


SECTION  III 


THE  NUMERICAL  DIFFERENCE  SCHEME 


In  this  section  of  the  report,  the  MICE  numerical  difference 
scheme  is  described.  The  starting  point  of  this  description  is  a 
presentation  of  the  set  of  partial  differential  equations  useJ  in  MICE 
for  the  flow  of  an  infinitely  conducting  fluid.  The  approximations 
used  in  the  difference  equations  are  then  given.  The  concept  of 
difference  operators  is  used  to  explain  wlmt  is  meant  by  time-step¬ 
splitting  and  the  MICE  difference  operator  is  then  described  in  suf¬ 
ficient  detail  so  that  the  interested  reader  can  understand  the  numerical 
scheme. 


Two  topics  not  discussed  in  this  report  arc  the  MICE  transport 
of  chemical  species  and  Lagrangian  coordinates.  Information  on  the 
Lagrangian  coordinates  is  contained  in  Reference  1  which  discusses  the 
idea  of  Lagrangian  coordinates  in  an  Eulorian  code  nnd  gives  a  r'irst 
order  scheme  for  advancing  the  quantities  in  time. 


DIFFERENTIAL  EQUATIONS 


3p 

at 


1^2  .  . 


(pv)  ri) 

(pvv)  -  v  (p  +  +(b  •  $)  t  -  PRez 
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(2) 


a! 

3t 

4x(vxl) 

(3) 

3 

3t 

(pi)  -  -  $  •  (pi  v)  -  V  ( 

$  .  v) 

(4) 

p 

-  (Y  -  1)  pi 

(5) 

where  p  is  the  mass  density  in  gm/cc,  V  is  the  tluid  velocity  in 
cm/sec,  P  is  the  fluid  pressure  in  dynes/cm*,  1  is  the  magnetic 
field  strength  in  >/4tT  gauss,  I  is  the  specific  internal  energy  in 
ergs/gm,  and  y  is  the  specific  heat  ratio. 

APPROXIMATIONS 

The  collection  of  nil  values  of  p,  pV,  B,  pi,  R( ,  Jn^|  at 
all  points  in  the  space  mesh  can  be  considered  to  be  the  components  of 
a  multi-dimensionnl  vector  U.  An  explicit  difference  approximation 
to  the  differential  equations  has  the  form 

UnM  "  U"  ♦  Lx(Un)  +  Ly(Un)  +  Lz(Un)  ,  (6) 

whore  L  ,  L  ,  L  are  operators  which  involvo  taking  differences  in 
x  y  z 

the  x,  y,  and  z  directions,  respectively,  and  the  superscripts  indicate 
the  time  level.  An  implicit  difference  approximation  has  the  more 
general  form 

Un+1  -  Un  +  Lx(Un,  Un+l)  +  Ly(Unf  Un+1)  +  Lz(Un,  Un+1)  (7) 

As  is  indicated,  the  difference  operations  involve  the  time  advanced 
(unknown)  quantities  as  well  as  known  quantities.  In  this  case,  the 
solution  of  the  difference  equations  for  Un+l  may  involve  solving  a 
very  large  set  of  coupled  algcbruic  equations.  This  is  usually  accom¬ 
plished  by  either  matrix  techniques  or  iterative  techniques.  The 
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approach  taken  in  MICE,  is  to  use  implicit  difference  operators,  and 

to  invoke  an  additional  approximation  known  as  timo-step-splitting.  In 

this  approximation,  the  difference  operators  L  ,  L  ,  and  L  are 

x  y  z 

applied  sequentially  instead  of  simultaneously.  That  is,  we  define 
intermediate  values  of  U,  which  we  si. all  call  Un+X,  and  un+x+y,  as 
follows : 


Step 

1: 

un+x 

=  Un  +  L^(Un,  Un+X) 

(8) 

Step 

2: 

un+x+y 

"  Un*X  ♦  Ly(Un+X,Un+X+y) 

(9) 

Step 

3: 

un+1 

»  Un*X+y  ♦  L2(Un+X+y,  Un+1) 

(10) 

The  error  generated  by  this  additional  approximation  of  timo-step- 
splitting  is  discussed  in  Appendix  B. 

Each  of  the  above  steps  involve  the  application  of  one  difference 
operator.  The  advantage,  from  a  coding  standpoint,  of  this  approach, 
is  that  a  single  generalized  "difference  operator"  subroutine  can  be 
written,  which  will  perform  any  of  the  three  stops  above.  The  advantage 
from  standpoint  of  core  requirements,  is  that  only  the  values  corresponding 
to  a  single  row  or  column  of  cells  need  be  in  main  core  at  a  given  time. 
Thus,  the  number  of  computational  cells  which  may  be  used  for  a  given 
problem  is  limited  not  by  main  core  size,  hut  by  the  size  of  the  F.CS 
bank  on  the  COC  6600,  or  LCM  on  the  7600.  One  further  advantage  is 
that  it  bocomes  relatively  easy  to  switch  between  doing  1-0,  2-0,  and 
3-0  problems. 

We  shall  now  present  the  differentia!  equations  to  which  the 
ono-dimonsional  difference  operator  must  correspond.  In  these  equations, 
the  are  transport  coefficients  which  urc  adjusted  to  improve  the 
numerical  stability  of  the  code,  Q  is  the  artificial  viscous  pressure, 


b  and  B  are  the  parallel  and  perpendicular  magnetic  field  components, 

u  and  V  are  the  parallel  and  perpendicular  velocity  components, 

P*  a  P  +  (B2+b2)/C2  is  the  jo-called  "Boris  mass",  and  m  =  0,  1.  or  2, 

Si 

depending  on  whether  the  difference  operator  is  being  applied  in  planar, 
cylindrical,  or  spherical  geometry. 


The  equations  are: 


^CP)  ■ 

/pu 

m  3x 

-  “>  ST  e) 

(ID 

(b) 

1  3  m  . 

*  -  u  — —  X 

m  3x  1 

(b  ’"IS1) 

(12) 

X 

“  *  -1  7x  *  (Bu  •  rt3  4  B)  +  b  ET  v 

(13) 

X  ' 

H(pu) 

1  3  ml 

-  "  “i  37  X  [P 

m‘  -  a,  o  u) 

x  ' 

- 

.  ~  CP+Q)  -  p/p* 

114) 

^  (  g  ,  for  vertical  direction 

whcrc  G  "  j  0,  otherwise 

Ji  fpVD“  “  'Tn  Sx  *  (pVu  ‘  a2  p  ^  V) 

+  p/p*  b  ^  D  (15) 


We  arc  considering  only  one  perpendicular  component  of  the  velocity 
and  B  field.  The  generalization  to  two  perpendicular  components  is 
straightforward. 
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3  1  3  m  (  3  T\ 

St  (pl>  "  -  -a  S  *  (plu  -  3x  pIj 

-  (p*«  -i  h  '  “  <16> 

x 

or, 

i  <E>  ■  -  ji  A  *”  (tE‘p*«lu  - s E)  • 0  G’ u 

(17) 

where  E  *  pi  +  ~-p  (u2+V2)  (18) 

The  code  may  be  run  using  difference  forms  of  either  equation 
(16)  or  equation  (17)  for  the  energy.  The  latter  choice  results  in 
better  energy  conservation,  but  equation  (16)  gives  more  reliable 
temperatures  in  regions  where  magnetic  forces  dominate.  The  factor 
p/p*  which  appears  in  front  of  the  magnetic  force  terms  has  the  effect 
of  limiting  the  Alfven  speed  to  the  value  C,  in  regions  where  the 
mass  density  is  very  small.  A  limiting  Alfven  speed  of  2xl07  cm/sec 
has  been  found  to  be  a  satisfactory  value. 

The  difference  operator  is  constructed  from  these  one 
dimensional  differential  equations.  Multi -dimensional  calculations  are 
performed  by  suitably  rotating  the  vector  components  in  the  routine 
which  fetches  the  cell  quantities  from  large  core  memory,  so  that  the 
operator  routine  always  "thinks"  it  is  going  in  the  x-direction. 

Lx  DIFFERENCE  OPERATION 

Another  considerable  simplification  in  the  equations  can  be 
made  by  further  splitting  this  set  of  equations  into  two  sets,  which 
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we  will  call  "L^  phase  1"  and  "L^  phase  2."  This  splitting  is  done 
on  the  b'.sis  of  the  distinction  between  transverse  Alfven  waves  and 
longitudinal  or  compressional  Alfven  waves.  Phase  1  calculates  the 
changes  in  the  perpendicular  components  of  B  due  to  the  propagation 
of  transverse  waves,  while  phase  2  calculates  the  changes  in  these 
components  of  B  due  to  propagation  of  compressional  waves.  Phase  2 
also  calculates  the  changes  in  all  other  quantities. 


The  phase  1  equations  are 


&<»  ■ 

bH 

V  , 

(19) 

and 

&  <pV) 

»  p/p* 

‘  b  f  B  . 
dx 

(20) 

For  the 

case  constant  b 

and  constant  p,  these  equations 

are  easily 

seen  to 

be  equivalent  to 

the  pair  of  2nd  order  equations 

a2 

bl 

0* 

a2 

(21) 

3 2 

32 

and 

3tr  a 

p* 

3**(v>  • 

(2 

These  are  wave  equations,  representing  waves  traveling  with  velocity 
+  1  '/yfp*  .  The  amplitudes,  Vo  and  B„,  are  related  by  the  expression 

Bo  ■  /P*  V0  .  (23) 

Each  application  of  the  difference  operator  results  in  the 
Phase  1  equations  being  solved  implicitly  for  the  (partially)  time 
advanced  values  of  B,  after  which  the  phase  2  equations  are  solved 
for  the  time  advanced  values  of  nil  of  the  cell  quantities. 
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L  PHASE  1  DIFFERENCING 

X 


The  difference  equations  used  in  phase  1  for  advancing  the 


values  of  B  are 


Brt+1  ■  Bn  + 
3  J 


r  .  vn+is  l 

xj  [  j+*s  .i -*4  J 


vn+l*  -  vn  ,  ♦ 


j+*S 


-  rr\  -  b"+1  . 

j+i  j 


(p >;♦»,  Axj 


where  the  integer  subscripts  refer  to  values  ut  cell  centers,  half- 
integer  subscripts  refer  to  values  at  cell  boundaries,  and  the  super¬ 
scripts  refer  to  the  time  level.  Cell  boundary  values  are  obtained  by 
interpolation  In  the  array  of  cell  center  values.  These  equations  are 
solved  for  the  quantities  B,  ,  which  are  then  used  as  Initial  ''ulues 
for  L  phase  2.  These  equations  dLffer  from  the  usual  Lux-Wendroff 

A 

two  step  form  in  that  the  right  hand  side  of  equation  (25)  contains 
the  tlme-udvanced  values  Bj+l .  To  obtain  the  solution,  equation (25) 
is  substituted  into  equation  (24),  giving  the  equation 


“1  sr  *  b  »';*! 


where 


CAt°2  hJ 

2Ax.  Ax,.  (P* 


j  +,4 


(At) 2  b "  b" 


2AXj  AXj  ^  (p  )ji>j4 


(29) 


and 


(j.  *1  -a.  -  y. 

3  .1  3 

At  bn 

W.  =  Bn  +  — T — 

j  j  Ax. 


(30) 


The  collection  of  equations  (26)  for  all  values  of  j,  togother 
with  the  boundary  conditions,  constitutes  a  "tridiagonal"  matrix,  the 
solution  to  which  is  discussed  in  Appendix  A.  The  boundary  conditions 
are  discussed  in  Appendix  C. 


Lx  PHASE  2  DIFFERENCING 

The  phase  2  equations  arc  equations  (11)  through  (18),  except 
that  the  second  term  on  the  right  hand  side  of  equation  (13)  is  omitted, 
its  effect  having  been  accounted  for  in  phase  1.  The  difference  equations 
for  phase  2  can  be  simplified  by  introducing  symbols  for  the  flux  of 
mass,  momentum,  etc.,  as  follows: 


!>*'  .  Pr 
i  •'  j 


n+ 1 
i 


(.31) 


.n+S  n+!; 


„n+i  n+i  n+ 1  n+i 

V  uj*i  ‘  pj 


(ov<-(pu)^v:;;h  [*;;;  v- .  v- 


(bf) 


-  bn ,  -  ^  aj 

J+*l  J+lt  2 


-  x  b 


x”'  ,  Ax.  . 

)+!i  J+l* 


,-u^n+S  _  <'a^j+ls  /u»+1  un+«  \ 

(bf  V-.  ■  •  (V-  ‘  bj  ) 


(32) 


(3.3) 


(34) 


(.34  J) 
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V  j+ls  j+H  j  '*•*4 


(as) 


Ax 


J±!t 
.HS  L 


n+i  n+i 

8  j*>  '  "j 


(3S) 


„n+?4  .  ,v  n  rtn+>4  n4»a 

Lj  +  's  (Yj^  ‘  n  pj^  V* 


r  * 

M +4  ,,n+i  ,,n+» 

‘  Ljtl  '  ‘J  J 


tj  ,  1  n+Jj 

+  i  u.  , 

s  j+JiJ  J+,s 


(36) 


with  E  ■  pi  ♦  i  p  (uVj1! 


n<l*  -  (pu)n+l»  il1+H  -  -l!±h. 

-  lPUJJ+»j  ‘j+lj  Av 


Ax 


J+,i 


n+i  n»i  iim  ,n+i. 

P.U>  1 J  +  i  *  °)  *j  I  {37) 


1  <’ 


Equations  (311  through  (37)  above  define  fluxes  at  the  cell  boundaries, 
Cell  cantered  fluxos  nre  defined  In  the  MICE  scheme  as: 


(Ouf)  ^ 

■  >;B)  ■ 

(38) 

(pVF)  ( 

pj  uj  vi  ’ 

(39) 

(Bf) 

-  n"u!f+l 

J  .1  » 

HO) 

(Pff)j 

-  p"+I  i!;  1.7  . 

(41) 

The  cell-boundary  fluxes  are  time-centered,  except  for  the 
additional  transport  terms  Involving  the  u's  and  differences  of  the 
time-advanced  quantities.  The  a's  havo  the  form 

ai  "  I  *  I  1 1  (Ax) 2 1  >  (42) 
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^jgggj|j|ggggggjjfigyi 


Since  these  terms  are  of  second  order  in  Ax,  we  see  that  the  cell¬ 
boundary  fluxes  are  time-centered,  to  first  order.  The  cell-centered 
fluxes,  however,  are  mixtures  of  time  n  and  time  n+i  quantities. 
We  therefore  refrain  from  assigning  a  time  level  to  the  symbols  for 
these  fluxes. 


We  shall  now  write  out  the  MICE  difference  equations,  using 
difference  operators  dofined  as: 


‘"’i  ‘  ‘ 


'J  j 


■ x)  Aj] 

Vs  AVS  L  J 


(43) 


(44) 


(D),A  5 


'j 


=  ,-L  Ta 
■  ^  . 


j+ij "  Vs] 


(lV.A2  [V>-Ai] 


(45) 

(46) 


The  difference  equations  are: 


n+1  _  _n 


Pj  •  At  (V»),  (Pf) 


n+S 


(47) 


.n+1  _  ,n 


n+l 


b”  -  At  u'jTl  (VO,  (bf )"+l*  -  At  (V-),  (bfO"^  (48) 


v^n+S 


-  Bj  -  At  (V*),  (Bf )n+**  , 


(49) 
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»  b* 

whore  a*  «  y  (y-1)  I  +  /*  is  the  effective  sound  speed  squared, 


(pn$  ■  (Pin 


[<V.  ‘  1)(l,nU  *  Qj*-,] 


(55) 


»"3  •  *"♦>,  •  (Bf>  ■ 


(56) 


u  ml 


(r>V)J^  ■  (oV)J4S  -  }a,(V-)ja  fnvfi 


‘*4,Wm  "U 


„n+ 1 


(57) 


The  solution  to  the  above  difference  equutians  is  obtuinod  ulgebvaieul ly 
substituting  equations  (54)  und  (51)  into  equation  (47) ,  obtaining  a 
tri-diagonal  set  of  equations  for  the  o  1 ,  which  can  then  bo  solved, 

The  quantities  p'|*^  and  are  then  calculated  by  Interpolation 

in  space  and  time.  Having  the  values  for  pl1+  1  ,  equation  (54)  Is 

then  used  to  obtain  (pu)'!*!4  ,  which,  together  with  i>?V  ,  then  gives 

)  +  J  J  *  5 

L 

the  valuos  for  , 


liquations  (32),  (40),  (41),  (55),  and  (56)  arc  combined 
algebraically  with  oquation  (50)  to  obtain  a  tri-diagonal  set  of 
equations  which  can  be  solved  numuricn 1 ly  for  the  quantities  u'|+ 1 , 

Equations  (55)  and  (56)  can  then  be  used  to  obtain  the  quantities 

(pl)”*jj|  ,  and  .  EquatJ  -ns  (48)  und  (49)  are  then  solved  to  find 

the  new  magnetic  field  comp v.nonts ,  after  which,  equation  (57),  and  then 
equation  (51),  may  bo  solved  for  the  now  perpendicular  velocity  components. 
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Then,  finally,  the  energy  equation  is  solved  (equation  (S2)  or  equation 
(iJ),  depending  on  which  form  is  desired.) 

When  solving  each  of  the  above  equations,  the  appropriate 
boundary  conditions  must  be  considered.  Reflective  boundaries  and 
poriodic  boundaries  are  relatively  straightforward,  and  are  not  dis¬ 
cussed  in  this  report.  Tranamittive  boundaries  are  more  troublesome, 
however,  and  the  MICE  treatment  of  such  boundaries  is  discussed  in 
Appendix  C. 
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APPENDIX  A 


THE  SOLUTION  OF  A  TR I -DIAGONAL 
SYSTEM  OF  LINEAR  EQUATIONS 


A  system  of  N  linear  equations 
can  be  put  in  the  form 


A.  II.  +  B.U.  +  C.U.  -  W. 
J  J  +  i  3  J  J  J-i  ) 


A  II  +  B  U 
12  11 


V’n  *  CNUN-t  ■  »N 


is  called  tri-diagonal  if  it 

2  i  j  <  N  -  1  (Al; 

(A2) 

(A3) 


whore  the  ll.'s  are  the  unknown  quantities,  and  A^  ,  B^  ,  ,  and 

W.  arc  known.  The  equations  may  lie  ordinary  scalar  equations,  or 
they  may  be  matrix  equations.  Where  equations  of  this  form  rise  in 
.  1  ICIi: .  equation  (Al)  represents  the  difference  equations  in  the  in  erior 
of  the  mesh,  while  (A2)  and  (A3)  are  the  boundary  conditions. 


The  solution  is  obtained  by  introducing  additional  quantities 

X  ,  Y.,  such  that 

.1  .1 

U.  =  X.U.  +  Y.  (A4) 

J  +  i  J  J  I 


Substituting  equation  (A4)  into  equation  (Al),  we  obtain 


(Aj  Xj  +  B.) 


+  C.U.  =  W.  -  A . Y . 
J  d-i  J  .1  } 


(A5) 
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Raising  the  index  j  by  one,  and  solving  for  Uj+1»  we  obtain 

U.  -  “  (a.  X.  +  B .  V 1  C.  U. 

3  +  i  \  j+i  j  +  i  J*i]  J*i  ] 

*  (Ai±  L  ♦  B ,  V1  “A*  \ 

\  j  +  l  j*i  J+t)  \  3+i  jvi  j*i ) 

Comparing  this  expression  with  equation  (A4),  we  make  the  identifications 

X,  "  "  (Au  x^  +  +  B<  +  V‘  cu  *  md  (A6) 

j  \  J+i  J+i  J+i)  J  +  i 

Y.  -  (A-.  +  B.  Yl  /W.  -  A.  Y.  \  . 

J  \  J*i  3*1  J*i J  ^  j+i  j*i  j+ij 

Comparing  equations  (A3)  and  (A4),  we  see  that 


'N- 


N’ 


,  ■  -  (®n)  '*  CN  '  and 

,  '  (bn)  ''  "n  ■ 


Equations  (A6)  then  give  the  X's  and  Y's  for  all  smaller  values  of  j  . 

Equations  (A41,  with  j  «  1,  is  then  combined  with  equation  (A2),  to 
give 

Ui  "  (AiXi  +  Bi)_1  (Wi  '  AiYi)  ’ 

The  Uj's  for  all  higher  values  of  j  can  then  be  found  from  (A4) . 
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APPENDIX  B 


ANALYSIS  OF  THE  ADDITIONAL  ERROR 
GENERATED  BY  TIME-STEP-SPLITTING 


For  simplicity,  let  us  consider  the  case  of  two-dimensional 
cartesian  geometry.  The  differential  equations  we  are  attempting  to 
solve  can  be  put  in  the  form 

lr  ■  (A  *  “  hr)  H  ■  (B1) 

where  the  components  of  the  vector  U  are  the  MUD  dynamical 
variahlos  P,  I,  Vx,  V  ,  ,  By,  and  A  and  B  are  matrix  functions 

of  the  components  of  U  .  Then, 


_r_ 

u  = 

/  A  - — 

*  B 

a  ^ 

1  _i_  [)  *  ( 

A  9-  ♦ 

bL' 

\2  u 

'it2 

\  >1X 

) 

}  ?t  -  \ 

9X 

*Yt 

1  ~ 

and 

dm 

u  * 

(a-4- 

+  B 

4  ] 

1 111  U  . 

(B21 

9tm 

\  9* 

3y ) 

A  difference  approximation  to  the  differential  equations  has  the 

form 

Un+1  =  n  (Un,Un+1)  (B3) 

where  SI  is  a  difference  operator,  and  the  superscripts  indicate  the 
time  level.  The  statement  that  0,  is  accurate  to  order  l  in  At  , 


T 


means,  mathematically,  that  if  each  term  of  the  right  hand  side  of  equation 
(B3)  is  expanded  in  a  Taylor  series  about  time  tn,  then  the  equation 
becomes 

,.n+i  ,.n  ^  /.8  ^  „3  \  „n  ^  At2  /.3  ^  \  '  ,,n 

U  -  U  r  At  *  B|y  j  U  *  —  (A37  +  f  U 


+  ...  +0  (At  1 )  .  ( 

’.\'e  say  then,  that  ■  I  +  At  +  B  j  +  +  ® 

Now,  let  us  suppose  that  we  are  given  difference  operators 
and  Ly,  which  correspond,  respectively,  to  the  differential  equations 

“St  -  ’  A  “5x  -  *  1 


TtH  "  B“5y 


A  non-time -split  difference  scheme  would  advance  the  values  of  U 
from  time  n  to  time  n+i  by  means  of  the  operation 


U  1  *  (L  +  L 

-  \  x  y 


)  («>"*>) 


In  MICE,  the  operatois  L  and  L  are  applied  sequentially, 

x  y 

rather  than  simultaneously,  however.  That  is, 


,,n+x  . 

U  *  L. 


(un,ur  X) 


Un+1  «  Ly  (un+X,Un+J)  -  Ly  (Lx(bn,Un+X)  ,  Un+‘) 


Now,  let  us  suppose  that  L  and  L  are  individually  accurate  to  second 

x  y 

order.  That  is. 
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2 


L 

x 


L 


y 


1  * At  A  I*  *  ~r  (A  k)  +  0 


1  +  At  B  fy  +  T 


(B  V) 


+  0 


(At3) 

(At3) 


The  "product"  of  these  two  operators,  keeping  terms  up  through  order 
At2  ,  is 


Tims,  wo  see  that  there  is  an  additional  error  introduced  by  the  time 
splitting.  This  error  is  of  second  order  in  At  ,  and  is  given  by 


I; 


(’ 


3 
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This  is  the  lowest  order  error  term  caused  by  the  time-splitting 
whm  the  order  of  application  of  I-x  and  l.  is  as  indicated  above.  If 
the  order  of  application  is  reversed,  then  E  -  E.  This  fact  makes 
possible  an  easy  check  on  the  effects  of  this  error  term.  To  this  end, 
two  fireball  rise  calculations  have  been  performed,  differing  from  each 
other  only  in  that  the  order  of  application  of  L  and  L  was  reversed. 

The  results  of  these  two  calculations  wore  essentially  identical,  indicating 
that  the  additional  error  generated  by  time-step-splitting  can  be  ignored- 


APPF.NDIX  C 


BOUNDARY  CONDITIONS 


The  typo  of  boundary  which  has  always  been  difficult  to  handle 
properly  is  the  out-flow,  or  "transmittive"  boundary.  Other  types  of 
boundaries  can  be  treated  in  relatively  straight-forward  ways. 
Therefore,  we  will  concentrate,  in  this  appendix,  on  the  MICE 
treatment  of  the  transmittive  boundary. 

At  points  which  are  outside  the  space  mesh,  we  nre  either 
not  interested  In  the  values  of  the  dynamical  variables,  or  we  decide 
not  to  attompt  to  compute  them,  on  the  basis  of  lack  of  computing 
capability.  In  any  ense,  since  the  size  of  computer  memories  is 
finite,  there  must  always  bo  an  "outor  boundary"  beyond  which  lie 
space  points  for  which  we  huve  no  dotailod  information.  Therefore, 
some  assumption  must  be  made  at  the  outer  boundary,  in  order  to  be 
able  to  compute  in  the  colls  adjacent  to  tho  boundary. 

In  the  MICE  code,  tho  following  assumptions  have  been  imposed 
at  the  outer  boundaries: 

(i)  There  can  be  no  "inflow”  of  mass. 

(ii)  The  time  rate  of  change  of  the  amplitude  of 
incoming  waves  is  zero. 

These  conditions  seem  to  allow  signals  to  propagate  out  of  the  mesh 
without  causing  excessive  disturbances  within  the  mesh. 

In  phase  1,  where  trar'.jrse  Alfven  waves  are  treated,  only 
condition  (ii)  above  is  applicable.  Let  us  consider  a  row  of  cells, 
numbered  from  1  to  JMAX,  with  J  increasing  from  left  to  right.  In 
coll  no.  JMAX,  we  want  the  tinr  rate  of  change  of  the  amplitude  of 
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the  left-going  wave  to  be  zero.  Now,  the  phase  1  differential  equations 
can  be  written 


and 


A.  /B,\»  B  — 

3t  \  V  x 

JL  /y  \  a 

at  (  i;  p* 


a 

37 


(Cl) 


(C2) 


Where  p*  ■  p  *  B*/C*  and  V|  „  B^  are  the  components  of  V  and 
B  perpendicular  to  the  x-diroction,  To  simplify  the  boundary  treatment, 
we  consider  Bx  and  p*to  he  constant,  Equations  (Cl)  and  (C2)  above 
ore  then  equivalent  to  the  pair  of  equations 

D 

£  (»,  *  \)  ■  k  (Bi  *  ^'vi)  ■  iC3> 

and  g 

±  (Bl  .  /p.v,)  .  4  (B,  -  ^Vt)  .  (C4) 


Equation  (C4)  states  that  constant  values  of  -  /p  travel  along 
B 

paths  with  X  -  t  *  constant,  while  (C5)  states  that  the 


B 


quantity  +  sfo*  V,  is  constant  along  X  +  . t  *  constant.  At 

1  Vp* 

the  right  hand  boundary,  for  the  case  B^>  0  ,  +  p*  Vj  is  the 
amplitude  of  the  incoming  wave.  Condition  (ii)  above  then  becomes 


=  0 


(C5) 
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Another  equation  if  needed  in  order  to  establish  the  rates  of  change  of 
Ba  and  VA  individually,  and  we  have  arbitrarily  chosen  to  use  the 
equation 


(C6) 


This  equation  is  implemented  in  the  followina  wav:  One  goes  hack 
into  the  mesh  a  distance  ^  ^  picks  up  tho  value  of  II  ^  from 

from  that  point,  and  uses  that  vnluo  for  the  new  Bj  at  the  right  hand 
boundary,  liquation  (C5)  then  determines  the  now  vnluo  of  VL 


This  treatment  producos  tho  wrong  value  for  the  amplitude 
of  the  outgoing  wave,  but  gots  the  desired  value  for  the  amplitude 
of  the  incoming  wave.  In  principle,  equation  (C4)  could  be  employed 
to  obtain  the  correct  value  for  the  amplitude  of  the  outgoing  wavo, 
but  experience  shows  that  attempts  to  assign  values  to  both  ineomin ; 
and  outgoing  waves  can  produce  values  at  the  boundary  which  nro  obviously 
non-physical.  Tho  method  we  have  chosen  is  loss  likoly  to  cause 
numerical  instabilities  at  the  boundary,  even  though  it  is  not  strictly 
mathematically  coi„-ect. 


The  abovo  is  a  description  of  the  right-hand  boundary  treatment 
for  Bx>0  .  Other  coses  are  treated  in  an  analogous  way. 

Tho  phase  2  equations, 


-|t-  (PD  +  (pV  *  0 

-Jt  (pyx)  +  (p  vx2  +  p‘)  *  0 


31 


4t  (I)  *  vx  4r  (I>  ♦  £  k  (v  -  0 


with 


P'  i  P  +  p/p4. 


»L 


and 


P  ■  (V-  1)  Pi  . 


can  bo  put  In  the  fos^i 


a 

"TT 

(' 

*  Vx 

••’nr  (• 

■  0 

(C7) 

;i 

17 

(n 

♦  Vnl  ♦ 

(V,  * 

,•1  h 
u  T.*x 

(o  +  vx)  ■  0 

(CH) 

3 

It 

(o 

*  V  * 

fVx  ■ 

in  t 

’  JX 

(u  -  Vx)  -  0 

(W) 

where 

(i  ■ 
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*  •! 
0* 

)"' 

and 

0  * 

f  JL 

J  PC 

, 

These  o«|U»t  1  onn  have  the  interpretation  timt  the  three  quanti  tl oh  I 
P  ♦  Vx,  and  a  -  Vx  ,  are  each  constant  ninny  character  1st l e  curves  whose 
slopes  arc 
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dx 

3t 


vx  4  c 


and 


dx 

nr  • 


respectively,  Again  confining  our  attention  to  the  right  hind  boundary, 
we  must  have  — ^  (o  •  Vx)  ■  0  ,  to  satisfy  condition  (11).  this 
causes  equation  (C9)  to  be  satisfied,  since  we  assume  that  the  amplitude 
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of  thv  incoming  wave  is  zero.  The  other  equations  which  we  use  at  the 
right  hand  boundary  are 


("St  *  vx-3»)  0  '  0  •  <C10> 

And 

(4r  *  vx  4r)  1  n  0  •  (Cll) 


If  k  and  I  satisfy  these  equations,  then  equation  (C7)  U 
satisfied.  Equation  (C8),  which  describes  outgoing  waves,  is  not 
satisfied,  but  this  is  apparent. Jv  not  as  important  as  meeting  the 
condition  of  no  incoming  waves. 

liquations  (CIO)  and  (Cll)  are  satisfied  by  moving  the  distance 
V^At  in  from  the  boundary,  picking  up  the  vnluea  of  p  and  I  from 
that  point,  and  using  these  as  the  new  boundary  values.  Dion,  the 
now  value  of  Vx  at  the  boundary  ia  chosen  to  satisfy  A(o  -  V  )  *  o  , 

with  the  restriction  that  V  in  the  boundary  coll  cannot  exceed  V 

X  X 

in  the  first  cell  inside  the  boundary. 

The  above  described  method  of  treating  transmitt ive  boundaries 
gives  satisfactory  results,  both  for  strong  shocks  and  weak  shocks 
passing  through  the  boundary. 
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